/*	
	This file is part of the Snoopy code.

    Snoopy code is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Snoopy code is distributed in the hope that it will be useful,
    but WITHOUT Ag.nyrtot WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Snoopy code.  If not, see <http://www.gnu.org/licenses/>.
*/



#include "common.h"
#ifdef MPI_SUPPORT
#include "transpose.h"
#endif
#include "debug.h"
#ifdef WITH_SHEAR


double time_shift(double t) {
	double tremap;
#ifdef TIME_DEPENDANT_SHEAR
	tremap = sin(param.omega_shear * t) / param.omega_shear;	// This is formally the angular displacement divded by param.shear= int dt S(t) /<S>
#else
	tremap = fmod(t + param.ly / (2.0 * param.shear * param.lx) , param.ly / (param.shear * param.lx)) - param.ly / (2.0 * param.shear * param.lx);
#endif
	return(tremap);
}

void remap(double complex qi[]) {
	int i, j, k;
	int nx, ny, nxtarget;
	
	DEBUG_START_FUNC;
	
	for( i = 0 ; i < g.nloc ; i++) {
		w1[i]=0.0;
	}
	
#ifdef MPI_SUPPORT
// We have to transpose the array to get the remap properly
	transpose_complex_XY(qi,qi);
	
	for( i = 0; i < g.nxtot; i++) {
		nx = fmod( i + (g.nxtot / 2) ,  g.nxtot ) - g.nxtot / 2 ;
		for( j = 0; j < g.nytot/g.nProc; j++) {
			ny = fmod( j + g.rankx * g.nytot / g.nProc + (g.nytot / 2) ,  g.nytot ) - g.nytot / 2 ;
			
			nxtarget = nx + ny;		// We have a negative shear, hence nx plus ny
			
			if( (nxtarget > -g.nxtot / 2) & (nxtarget < g.nxtot/2)) {
			
				if ( nxtarget <0 ) nxtarget = nxtarget + g.nxtot;
			
				for( k = 0; k < g.nztot; k++) {
					w1[k + g.nztot * nxtarget + g.nztot * g.nxtot * j] = qi[ k + i * g.nztot + j * g.nztot * g.nxtot];
				}
			}
		}
	}
	
	// transpose back
	transpose_complex_YX(w1,w1);

#else
	for( i = 0; i < g.nxtot; i++) {
		nx = fmod( i + (g.nxtot / 2) ,  g.nxtot ) - g.nxtot / 2 ;
		for( j = 0; j < g.nytot; j++) {
			ny = fmod( j + (g.nytot / 2) ,  g.nytot ) - g.nytot / 2 ;
			
			nxtarget = nx + ny;		// We have a negative shear, hence nx plus ny
			
			if( (nxtarget > -g.nxtot / 2) & (nxtarget < g.nxtot/2)) {
			
				if ( nxtarget <0 ) nxtarget = nxtarget + g.nxtot;
			
				for( k = 0; k < g.nztot; k++) {
					w1[k + g.nztot * j + g.nztot * g.nytot * nxtarget] = qi[ IDX3D ];
				
				}
			}
		}
	}
#endif

	for( i = 0 ; i < g.nloc ; i++) {
		qi[i] = w1[i] * g.mask[i];
	}

	DEBUG_END_FUNC;
	
	return;
}

void kvolve(const double tremap) {
	int i, j, k;
#ifdef _OPENMP
	#pragma omp parallel for private(i,j,k) schedule(static)	
#endif	
	for( i = 0; i < g.nxloc; i++) {
		for( j = 0; j < g.nyloc; j++) {
			for( k = 0; k < g.nzloc; k++) {
				g.kxt[ IDX3D ] = g.kx[ IDX3D ] + tremap * param.shear * g.ky[ IDX3D ];
			
				g.k2t[ IDX3D ] = g.kxt[IDX3D] * g.kxt[IDX3D] +
						       g.ky[IDX3D] * g.ky[IDX3D]+
							   g.kz[IDX3D] * g.kz[IDX3D];
							  
				if ( g.k2t[IDX3D] == 0.0 ) g.ik2t[IDX3D] = 1.0;
				else	g.ik2t[IDX3D] = 1.0 / g.k2t[IDX3D];
			}
		}
	}

	return;
}
#endif
